Preliminary Transcriptome Analysis of Long Noncoding RNA in Hypothalamic-Pituitary-Mammary Gland Axis of Dairy Cows under Heat Stress

Heat stress (HS) is directly correlated to mammary gland dysfunction in dairy cows, especially in summer. The hypothalamic−pituitary−mammary gland axis (HPM axis) plays an important role in the regulation of stress response and lactation physiology in heat−stressed dairy cows. The aim of this study was to explore the lncRNA profile, and the competitive endogenous RNA (ceRNA) regulatory network in hypothalamus, pituitary, and mammary gland tissues of heat−stressed and normal dairy cows. We performed RNA sequencing (RNA−seq) to identify differentially expressed (DE) lncRNAs, and the ceRNA regulatory network was established in HPM−axis−related tissues. Our results showed that a total of 13, 702 and 202 DE lncRNAs were identified in hypothalamus, pituitary, and mammary glands, respectively. Of lncRNAs, 8, 209 and 45 were up−regulated, and 5, 493 and 157 lncRNAs were down−regulated. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses indicated that DE lncRNAs target genes that might play a role in hormone synthesis, secretion and action, apoptosis, mitogen−activated protein kinase (MAPK), AMP−activated protein kinase (AMPK), and mechanistic target of rapamycin (mTOR) signaling pathway. Moreover, the ceRNA regulatory network associated with the MAPK signaling pathway in HPM−axis−related tissues contains 3286 lncRNA–mRNA pairs. Furthermore, the ceRNA regulatory network associated with apoptosis, prolactin, AMPK, and mTOR signaling pathway in the mammary gland contains 772 lncRNA–mRNA pairs. Thus, some lncRNAs may be involved in the regulation of stress response and the physiological process of lactation. The changes in lncRNA expression profiles and ceRNAs (lncRNA–miRNA–mRNA) in HPM−axis−related tissues are the key to affect the stress response and lactation physiology of dairy cows under HS, which provide a theoretical basis for the molecular mechanism in the stress response of HPM−axis−related tissues in dairy cows under HS.


Introduction
A series of stress reactions of dairy cows caused by high−temperature environments have been widely studied. The cumulative research shows that heat stress causes the reduction of feed intake, lactation performance, and reproductive performance of dairy cows, as well as immune dysfunction, which directly and/or indirectly affects the health of dairy cows [1]. As a major production index of dairy cows, lactation performance is easily affected by heat stress in summer [2]. Presently, cows cope with heat stress in summer mainly through physical cooling and dietary nutrition changes [3,4]. However, the hypothalamic-pituitary-endocrine axis regulates the homeostasis imbalance caused by environmental thermal stimulation [5]. Among them, hypothalamic-pituitary-mammary gland axis (HPM axis) plays an essential role in stress response and physiological regulation of lactation [6]. For example, prolactin (PRL) and growth hormone (GH) secreted by the pituitary gland are closely related to the development and lactation of mammary glands. Accordingly, it is particularly important to analyze the molecular mechanism of heat stress regulating the HPM axis in dairy cows for improving their lactation performance.
Long noncoding RNA (lncRNA) is a type of RNA that does not have the ability to encode proteins and is no more than 200 nucleotides in length. As an important part of gene transcription, lncRNA can affect the expression of functional genes through transcriptional regulation and posttranscriptional regulation [7]. Recently, lncRNA has been considered as having an essential role in mammary gland development and lactation [8]. Zheng et al. profiled the transcriptome of the mammary gland using RNA sequencing (RNA−seq) technology at peak and late lactation stages of Holstein cows, and 12 differentially expressed (DE) lncRNA potentially played important roles in bovine lactation [9]. The lncRNA mammary proliferation and fatty acid synthesis−associated transcript (MP-FAST) is highly expressed in the Holstein cow mammary gland during the middle lactation period compared to the dry period, and lncRNA MPFAST regulates the expression of the genes in the phosphoinositide 3−kinase/protein kinase B (PI3K−AKT) signaling pathway through sponging miR−103 and promotes the proliferation and synthesis of fatty acids of BMECs [10]. In the past few years, accumulating evidence has indicated lncRNAs are widely involved in regulating heat stress response in animals [11,12]. Researchers found that 24,795 novel and 3763 known lncRNAs were expressed in bovine mammary gland by using RNA−seq in heat−stressed and nonheat−stressed cows [13]. However, the comprehensive analysis of the differential expression of lncRNA in the HPM axis of dairy cows under heat stress has not been reported.
Therefore, this study used RNA−seq to analyze the differential expression profile of lncRNA in the HPM axis under heat stress and determined and clarified the role of lncRNA in regulating heat−stressed dairy cows. This result provides a theoretical basis for studying the molecular mechanism of the HPM axis.

Animal Ethics Statement
All experimental procedures were carried out following the Institutional Animal Care and Use Committee of Nanjing Agricultural University, China (SYXK 2011−0036).

Animals and Sample Collection
There was continuous monitoring of temperature and humidity index (THI) of the cowshed, rectal temperature, and respiratory rate of Holstein cows under the same management level in summer (August, heat−stress group (HS), n = 20) and winter (December, nonheat−stress group (NHS), n = 20) in Jielong Animal Husbandry Co., Ltd. (Huai'an, Jiangsu, China). Information about parity, lactation days, and milk yield of the cows in the NHS group and HS group is shown in previous articles [14]. In brief, the temperature (T, • C) and relative humidity (RH, %) of the cowshed were measured at 10:00 and 20:00 every day for a month. Electronic thermometers at the six fixed points in the cowshed were read, respectively. THI criteria were derived from the NRC (1971), and THI was calculated according to THI = (1.8 × T + 32) -(0.55 -0.55 × RH) × (1.8 × T − 26), where T is the temperature ( • C) and RH is the relative humidity (%) [15]. Rectal temperature and respiratory rate were measured by animal thermometer and stopwatch, respectively, at 10:00 and 20:00 every day of the experiment.

Detection of Heat Shock Proteins (HSPs) and Antioxidant Enzyme Activity
As a marker protein of heat stress, we detected the mRNA expression of HSPs in HPM−axis−related tissues, including heat shock protein 27 (HSP27), heat shock protein 70 (HSP70), and heat shock protein 90 (HSP90). Commercial kits were used to measure concentrations of malondialdehyde (MDA, Jiancheng Biological Project, cat: A003−1−2, Nanajing, China), superoxide dismutase (SOD, Jiancheng Biological Project, cat: A001−3−1, Nanjing, China), and glutathione peroxidase (GPx, Beyotime, cat: S0058, Shanghai, China) in HPM axis tissues according to the manufacturer's instructions. The activity of these antioxidant enzymes was calculated based on the absorbance. The experiment was repeated three times. Statistical analysis was performed using SPSS version 20 (SPSS Inc., Chicago, IL, USA). The significance of differences was determined using Student's t−test, significance was set at p < 0.05, and all data were expressed as the mean ± SEM.

Total RNA Extraction and Sequencing
We selected three tissue samples from each group for lncRNA sequencing. Total RNA was extracted using a Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. After total RNA was extracted, the rRNAs were depleted for the construction of a total RNA sequencing library. The digested products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using an Illumina Nova−Seq 6000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Long Noncoding RNA Identification and Differential Expression Analysis
To obtain high−quality clean reads, raw reads were further filtered by fastp (v0.18.0). Clean reads were compared to the cow reference genomes with HiSAT2 (v2.1.0), and the transcripts were reconstructed with the software StringTie (v1.3.4, Computer Biology Center of Johns Hopkins University, Baltimore, MD, USA) to obtain known transcripts and new transcripts. The CNCI (v2) and CPC (v0.9−r2) (http://cpc.cbi.pku.edu.cn/ (accessed on 15 September 2021)) were used to assess the protein−coding potential of novel transcripts by default parameters. The intersection of both nonprotein−coding potential results was chosen as lncRNAs. lncRNA differential expression analysis was performed by DESeq2 software between the two different groups. To identify differentially expressed genes (DEGs) across samples or groups, the edgeR package (http://www.r$-$project.org/ (accessed on 15 September 2021)) was used. We identified differentially expressed transcripts with |log2fold change(FC)| > 1 and false discovery rate (FDR) < 0.05.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) Validation
Total RNA was extracted with Trizol reagent (Invitrogen, Carlsbad, CA, USA). The density of samples' RNA was quantified spectrophotometrically at 260/280 nm. The RNA was reverse transcribed using Prime ScriptTM RT Master Mix (TaKaRa, cat: RR036A, Otsu, Japan) according to the manufacturer's protocols. The primers were designed with the NCBI Primer BLAST tool (https://www.ncbi.nlm.nih.gov/tools/primerblast (accessed on 1 August 2022)). The sequences of the primers of mRNAs and lncRNAs are listed in Table 1. All quantitative reverse transcriptase PCR (qRT−PCR) runs were performed on a StepOne Plus Real−Time PCR System (Applied Biosystems, Foster City, CA, USA) using ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd., cat: Q711−02, Nanjing, China) in a reaction volume of 20 µL. The relative quantifications of mRNAs and lncRNAs were normalized to β−actin with the 2 −∆∆Ct method.

Function Enrichment Analysis
We selected the differential expression analyses of the lncRNAs. Gene ontology (GO) enrichment analysis and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed to investigate the biological function of the target genes of DE lncRNAs using DAVID (http://david.abcc.ncifcrf.gov/ (accessed on 10 October 2021)).

lncRNA-miRNA-mRNA Network Association Analysis
Based on our previous transcriptome sequencing results and small RNA sequencing results, we established ceRNA regulatory networks for DE lncRNAs. The ceRNA network was constructed based on ceRNA theory as follows: (1) the expression correlation between mRNA-miRNA or lncRNA-miRNA was evaluated using the Spearman rank correlation coefficient (SCC). Pairs with SCC <−0.7 were selected as negatively coexpressed lncRNA-miRNA pairs or mRNA-miRNA pairs, and both mRNA and lncRNA were miRNA target genes, and all RNAs were differentially expressed. (2) The expression correlation between lncRNA-mRNA was evaluated using the Pearson correlation coefficient (PCC). Pairs with PCC >0.9 were selected as coexpressed lncRNA-mRNA pairs, and both mRNA and lncRNA in these pairs were targeted and negatively coexpressed with a common miRNA. (3) A hypergeometric cumulative distribution function test was used to test whether the common miRNA sponges between the two genes were significant. As a result, only the gene pairs with a p−value less than 0.05 were selected.

Comparison of Environmental Temperature and Humidity Index (THI) and Cow Information
Researchers have shown that THI < 72 has no effect on dairy cows. It could cause slight heat stress in dairy cows between 72 and 78, moderate heat stress in dairy cows between 78 and 89, and severe heat stress in dairy cows above 90 [16]. In our study, the THI value of the HS group was significantly higher than that of the NHS group ( Figure 1A, p < 0.001). Meanwhile, RT and RR of cows in the HS group were significantly higher than those in the NHS group ( Figure 1B,C, p < 0.001). These results showed that the cows in the HS group were in a state of heat stress relative to the NHS group.

Effect of Heat Stress on Heat Shock Proteins and Antioxidant Enzymes in HPM-Axis-Related Tissues of Dairy Cows
To evaluate the effect of heat stress on HPM−axis−related tissues of dairy cows, we detected the mRNA expression of heat shock protein. We found that the expression of HSP27 in HPM−axis−related tissues had no significant effect between NHS and HS groups ( Figure 2A, p > 0.05). The expression of HSP70 was up−regulated ( Figure 2B, p < 0.01) in HPM axis tissues of heat−stressed dairy cows. Moreover, the expression of HSP90 was significantly increased in pituitary and mammary glands (p < 0.01), but it had no significant effect in the hypothalamus (p > 0.05) under heat stress conditions ( Figure 2C). Subsequently, we directly detected the content of antioxidant enzymes in HPM axis tissues. The results show that MDA content was significantly increased in HPM axis tissues between NHS and HS groups ( Figure 2D, p < 0.01). SOD activity significantly decreased ( Figure 2E, p < 0.01) in the hypothalamus and mammary gland, and GPx activity significantly decreased ( Figure 2F, p < 0.01) in pituitary and mammary glands under heat−stress conditions.

Identification of lncRNAs in HPM-Axis-Related Tissues of Dairy Cows under Heat Stress
In the current study, 18 cDNA libraries in the HPM−axis−related tissues were constructed to identify the lncRNAs expressed in the NHS and HS groups. After the raw reads were subjected to quality control and filtering, high−quality filtered effective reads were obtained from the 18 libraries. The filtered reads were subsequently compared with the cow reference genome, and 81.08-85.63% reads were successfully mapped to the reference genome. Therefore, the DEG analysis based on the genome was considered reliable. A total of 2163 lncRNAs were detected in the NHS and HS, and 1538 novel lncRNAs were screened on the basis of the data from the CNCI and CPC2 programs ( Figure 3A). We statistically analyzed the numbers of lncRNAs on the basis of their specific positions and found that intergenic lncRNAs accounted for 70.22% of all the lncRNAs ( Figure 3B).

Differential Expression Analysis of lncRNAs
As shown in Figure 4A, we compared the expression abundance of all transcripts in 18 libraries. In general, the relative expression levels were lower in the NHS group than in the HS group, and the trend of down−regulation was obvious. The overall expression levels of HPM−axis−related tissues were also compared in the NHS and the HS groups ( Figure 4B), and the expression levels in the HS group were slightly lower than those in the NHS group. Further, a total of 13, 702 and 202 DE lncRNAs were identified in the hypothalamus, pituitary, and mammary gland, respectively. Of lncRNAs, 8, 209 and 45 were up−regulated, and 5, 493 and 157 lncRNAs were down−regulated ( Figure 4C-F).

Validation of DE lncRNAs with qRT−PCR
We randomly selected three lncRNAs from HPM−axis−related tissues and detected their expression by qRT−PCR. The trends in the lncRNA expression changes were similar to those identified via RNA−seq, although some quantitative differences were found between the qRT−PCR and RNA−seq results ( Figure 5). These findings indicated that the RNA−seq results were reliable and could be used for bioinformatic analysis.

Function Enrichment Analysis for Differentially Expressed lncRNAs
In order to determine the biological processes associated with the target genes of DE lncRNAs in HPM axis tissues, we performed GO and KEGG pathway enrichment analysis. In the hypothalamus, DE lncRNA target genes were significantly enriched in 504 GO pathways (p < 0.01), among which 357 GO pathways were associated with biological processes, 88 with cellular components, and 59 with molecular functions ( Figure 6A, Table S1). The functional annotation by GO indicated that the predicted target genes in the hypothalamus were related to regulation of neurotransmitter levels, regulation of peptide hormone secretion, hormone transport, as well as multicellular organismal response to stress ( Figure 6B). In the pituitary, a total of 719, 173, and 162 GO terms (p < 0.01) were significantly enriched in biological processes, cellular components, and molecular functions, respectively ( Figure 6C, Table S2). The functional annotation by GO indicated that the DE lncRNA target genes were related to cellular response to hormone stimulus, regulation of corticosteroid hormone secretion, multicellular organismal response to stress, as well as positive regulation of translation in response to stress ( Figure 6D). Moreover, DE lncRNA target genes were significantly enriched in 1432 GO terms (p < 0.01), among which 263 GO pathways were associated with biological processes, 108 with cellular components, and 115 with molecular functions in the mammary gland ( Figure 6E, Table S3). The functional annotation by GO indicated that the DE lncRNA target genes were related to hormone−mediated signaling pathway, hormone receptor binding, and regulation of cellular response to stress ( Figure 6F). Furthermore, target genes of DE lncRNAs participate in energy metabolism, lipid metabolism, amino acid metabolism, signal transduction, cell growth and death, and environmental adaptation in the HPM−related tissues according to KEGG pathway enrichment analysis ( Figure 7A-C). A total of 35, 104, and 52 KEGG pathways were significantly enriched in the HPM−related tissues, respectively (p < 0.05, Tables S4-S6). We found that the DE lncRNA target genes in the hypothalamus were enriched in glutamatergic synapse, cyclic adenosine monophosphate (cAMP) signaling pathway, mitogen−activated protein kinase (MAPK) signaling pathway, and autophagy ( Figure 7D). In the pituitary, the DE lncRNA target genes were enriched in hormone−related pathways, such as growth hormone synthesis, secretion and action, thyroid hormone synthesis, cortisol synthesis and secretion, and prolactin signaling pathway ( Figure 7E). The KEGG enrichment analysis of DE lncRNAs in the mammary gland showed that the meaningful pathways were related to apoptosis, such as the MAPK signaling pathway. Simultaneously, DE lncRNAs are also enriched in lactation−related pathways, such as mechanistic target of rapamycin (mTOR) signaling pathway, AMP−activated protein kinase (AMPK) signaling pathway, and prolactin signaling pathway ( Figure 7F). From these analyses, we speculated that lncRNAs may participate in the stress response and physiological lactation process of dairy cows by regulating the expression of target genes.

Discussion
In a previous study, we measured the serum endocrine hormones and biochemical indexes of cows in an NHS group and HS group. Heat stress significantly increased the levels of adrenocorticotropic hormone (ACTH) and cortisol (COR) in serum of dairy cows, but decreased the levels of triiodothyronine (T3) and thyroid hormone (T4). These results showed that the stimulation of the thermal environment caused the changes in the hypothalamic-pituitary-endocrine axis in dairy cows [14]. Further, in the current study, we used the transcriptome to analyze the differential expression profile of lncRNA in the HPM axis of heat−stressed dairy cows, and established the regulatory network of lncRNA-miRNA-mRNA in the hypothalamus, pituitary, and mammary gland. Due to limited funding, we could only randomly select the hypothalamus, pituitary, and mammary gland tissues of three cows from each group for transcriptome sequencing.
Heat shock proteins (HSPs) are polypeptides with a highly conserved structure, and are also called heat stress proteins since almost all organisms express heat shock protein when they are subject to extra stresses to a certain extent [17]. In order to determine the occurrence of heat stress in dairy cows, we detected the mRNA expression of HSPs in HPM−axis−related tissues. Our experiment found that HSP70 was significantly up−regulated in hypothalamus, pituitary, and mammary gland tissues of heat−stressed dairy cows. Meanwhile, the expression of HSP90 was significantly increased in pituitary and mammary glands. However, thermal stimulation did not cause a change in HSP27 in HPM−axis−related tissues. HSP70 and HSP90, as ATP−dependent molecular chaperones, play a role in helping the folding of newly synthesized peptides, the assembly of multiprotein complexes, and the transport of proteins across cell membranes [18]. However, HSP27 has multiple phosphorylation sites and plays its role through protein phosphorylation [19]. Therefore, these results may be related to the molecular structure of HSPs. In addition, previous studies have shown that heat stress can cause antioxidant dysfunction in dairy cows, which is mainly manifested in the reduction of antioxidant enzyme activities, such as SOD and glutathione peroxidase (GSH−Px), and the increase in lipid peroxide MDA [20]. Some researchers mitigated the heat stress response of dairy cows by increasing the activity of antioxidant enzymes to improve the lactation performance [21,22]. In the present experiment, we found that MDA content was significantly increased in HPM axis tissues between NHS and HS groups. Meanwhile, SOD activity significantly decreased in the hypothalamus and mammary gland, and GPx activity significantly decreased in pituitary and mammary glands under heat stress conditions. These results demonstrated that the antioxidant capacity of HPM−axis−related tissues may be reduced, which indirectly proved that cows were in a heat−stress state. Therefore, the tissue samples we collected can be used for subsequent sequencing studies.
The mechanisms of noncoding RNAs (ncRNAs) in lactation and stress response have been studied, and they play an important role in regulating RNA transcription, protein translation, and protein interaction. Studies have confirmed that lncRNA plays an important role in lactation and heat stress response. In this study, 81.08-85.63% reads were successfully mapped to the reference genome through high−throughput sequencing. Further, a total of 13, 702, and 202 DE lncRNAs were identified from the hypothalamus, pituitary, and mammary gland, respectively. We selected 12 DE lncRNAs from HPM−axis−related tissues for validation. qRT−PCR results were highly consistent with those of RNA−seq.
These results indicate that lncRNA may play a key role in the HPM axis of dairy cows after thermal environment stimulation.
Differential lncRNA expression is often accompanied by changes in cellular function, and lncRNAs can interact with coding genes in both cis and trans manners to fulfill their functions [7]. We performed GO analysis and identified the biological processes related to stress response processes of DE lncRNA target genes in three tissues. The main functions include cellular response to stress, regulation of cellular response to stress, stress−activated MAPK cascade, regulation of response to stress, and stress−activated protein kinase signaling cascade. Subsequently, KEGG pathway enrichment analysis showed that target genes of DE lncRNAs were closely related to stress response and physiological lactation processes. For instance, the target genes of DE lncRNAs in the hypothalamus were highly enriched in the oxytocin signaling pathway, gonadotropin−releasing hormone (GnRH) secretion, and insulin secretion. In the pituitary, the target genes of DE lncRNAs were also mainly enriched in hormone−related pathways, such as cortisol synthesis and secretion, thyroid hormone synthesis, growth hormone synthesis, secretion and action, and prolactin signaling pathway. COR, as the most sensitive hormone in the stress response, increases significantly after heat stimulation [23,24]. Researchers also found that the effects of short−term and long−term environmental heat on COR were obvious, initially increased due to acute stressors, and plasma levels decreased after long−term exposure to stressors [25]. Moreover, thyroid hormone plays a key role in regulating body temperature, energy intake, and metabolism [26]. The stimulation of the thermal environment can significantly reduce the secretion of thyroid hormone to regulate the thermal homeostasis of dairy cows [27]. Universally, prolactin (PRL) and growth hormone (GH) are essential in maintaining the physiological process of lactation. The study has shown that PRL and GH in hormone−response pathways are involved in energy metabolism during thermoregulation processes in cattle, and single nucleotide polymorphisms (SNPs) within genes of the PRL and GH pathways are predictors of reproductive phenotype in heat−stressed Holstein cows [28]. Thus, these results indicate that lncRNA may be involved in regulating the expression of hormone−related target genes in response to thermal stimulation. Simultaneously, the meaningful pathways of DE lncRNA target genes in the mammary gland were related to apoptosis, mTOR signaling pathway, and AMPK signaling pathway. Chen et al. found that heat stress inhibits the proliferation of BMECs, aggravates cellular oxidative stress, causes mitochondrial dysfunction, and promotes BMEC apoptosis [29]. Meanwhile, sirtuin 3 (Sirt3) can protect BMECs from heat stress damage through the AMPK signaling pathway [30]. These studies support our findings that lncRNAs may play an important role in heat−stress−induced cell damage. In addition, the AMPK−mTOR pathway is closely related to cell energy and material metabolism [31]. Some studies have shown that the AMPK−mTOR pathway is involved in regulating the synthesis of casein in mammary cells [31,32]. Liang et al. analyzed and predicted the lncRNAs related to milk fat metabolism, and found that the target genes of these lncRNAs were also significantly enriched in MAPK, AMPK, and mTOR signaling pathways [33]. Interestingly, we found that DE lncRNAs were significantly enriched in the MAPK signaling pathway in the hypothalamus, pituitary, and mammary gland in the current study. The MAPK family is a group of serine/threonine kinases, which can be activated by cytokines, neurotransmitters, hormones, and cell stress [34]. Accumulating evidence has indicated that thermal stimulation induces MAPK activation, thereby affecting cell proliferation and apoptosis [35,36]. Consequently, we speculate that the differential expression of lncRNAs in the HPM axis may play an important role in regulating the lactation process of heat−stressed dairy cows. These possible associations should be further explored in future study.
MicroRNA (miRNA), as a highly conserved small molecule ncRNA, can regulate various physiological and pathological processes [37]. Many studies have shown that lncRNAs have miRNA−response elements (MREs), which can act as competitive endogenous RNA (ceRNA) to regulate miRNAs and mRNAs, thereby exerting their biological functions [38,39]. For example, lncRNA AK017368 promotes proliferation and inhibits differentiation of myoblast cells by attenuating the function of miR−30c [40]. Liu et al. found that lncRNA−MEG3 is highly expressed in bovine skeletal muscle tissue, and lncRNA−MEG3 can promote bovine skeletal muscle differentiation by interacting with miRNA−135 and myocyte enhancer factor 2C (MEF2C) [41]. Additionally, CTTN−IT1 is a novel lncRNA that is a ceRNA of miR−29a and can promote skeletal muscle satellite cell proliferation and differentiation by restoring the expression of Yes−associated protein 1 (YAP1) when it is inhibited by miR−29a in Hu sheep [42]. In our present study, some ceRNA regulatory networks were established in the HPM axis of dairy cows under heat stress, such as MSTRG.  [44,45]. The Janus kinase-signal transducer and activator of transcription (JAK−STAT) signaling pathway and suppressors of cytokine signaling (SOCS) family genes play a crucial role in controlling cytokine signals in the mammary gland and thus mammary gland development [46]. Accordingly, we speculated that cows may experience changes in the main lactating pathway to respond to the thermal environment by regulating the expression of lncRNA and ceRNA (lncRNA-miRNA-mRNA) networks in the HPM axis. The specific mechanism should be studied in future research.
In conclusion, the current study indicated that the expression profiles of lncRNAs in HPM−axis−related tissues of heat−stressed and nonheat−stressed cows were different. Moreover, cows may experience changes in the main lactating pathway to respond to the thermal environment by regulating the ceRNA (lncRNA-miRNA-mRNA) networks in the HPM axis. Consequently, these results provide a theoretical basis for the molecular mechanism of the stress response and lactation physiology of HPM−axis−related tissues under HS. Funding: This study was supported by the Technical System of Modern Agriculture (Dairy Cow) Industry in Jiangsu Province (JATS (2022) 469). The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.
Institutional Review Board Statement: All experimental procedures in present study were approved by the Institutional Animal care and Use Committee of Nanjing Agricultural University. The procedures used in the present study were in compliance with the laws and regulations of China (based on the Institutional Animal care and Use Committee of Nanjing Agricultural University, China (SYXK2011−0036)), and the internationally accepted principles and guidelines for the care and use of agricultural, laboratory, or experimental animals. All methods are reported in accordance with ARRIVE guidelines for the reporting of animal experiments.

Informed Consent Statement: Not applicable.
Data Availability Statement: All raw transcriptome data reported in current article have been deposited at the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/ sra/PRJNA885909 (accessed on 30 September 2022)).